from __future__ import division
from sympy import *
x, y = symbols('xy')
h = symbols('h', integer=True)
import numpy as np
from numpy.numarray import Float64
#k, m, n = symbols('kmn', integer=True)
#f, g, h = map(Function, 'fgh')
#voi- variable of integration
grad_phi = [['-1/h', '0'], #0
            ['0', '-1/h'], #1
            ['1/h','-1/h'],#2
            ['1/h', '0'],  #3
            ['0', '1/h'],  #4
            ['-1/h', '1/h'],#5
            
            ['-1/h', '0'], 
            ['0', '-sqrt(3)/h'],
            ['1/h','-sqrt(3)/h'],
            ['1/h', '0'],
            ['0', 'sqrt(3)/h'],
            ['-1/h', 'sqrt(3)/h']]
#inner limits, outer limits
limits = [['(y, 0, x)','(x, 0, h)'], #0
          ['(y, x, h)','(x, 0, h)'],
          ['(y, 0, x+h)','(x, -h, 0)'],
          ['(y, x, 0)','(x, -h, 0)'],
          ['(y, -h, x)', '(x, -h, 0)'],
          ['(y, x-h ,0)', '(x, 0, h)'],
          
          ['(y, 0, x/sqrt(3))','(x, 0, h)'], #0
          ['(y, x/sqrt(3), h/sqrt(3))','(x, 0, h)'],
          ['(y, 0, (x+h)/sqrt(3))','(x, -h, 0)'],
          ['(y, x/sqrt(3), 0)','(x, -h, 0)'],
          ['(y, -h/sqrt(3), x/sqrt(3))', '(x, -h, 0)'],
          ['(y, (x-h)/sqrt(3) ,0)', '(x, 0, h)'],
          ]

triangle_accordance = [
                        [(0, 2), (5, 3)], 
                        [(0, 4), (1, 3)],
                        [(1, 5), (2, 4)],
                        ]



for ind in xrange(3):
    t1, t2 = triangle_accordance[ind][0]
    t3, t4 = triangle_accordance[ind][1]
    triangle_accordance.append([(t2, t1), (t4, t3)])
    
def get_according_triangle(i, triangles):
    result = []
    for j in range(2):    
        t1, t2 = triangle_accordance[i][j] 
        if t1 in triangles:
            result.append((t1, t2))
        elif t1 + 6 in triangles:
            result.append((t1+6, t2+6))         
    return result

def get_adjoining_triangles(i, triangles):
    result = []
    loop = range(6)
    for j in range(i-1,i+1):
        if loop[j] in triangles:
            result.append(loop[j])
        elif (loop[j] + 6) in triangles:
            result.append(loop[j] + 6)
    return result

#def get_opposite_triangle(i, n_triangle):
#    if n_triangle < 6:
#        t11, t12 = triangle_accordance[i][0]
#        t21, t22 = triangle_accordance[i][1]
#        if t11 == n_triangle:
#            return t12
#        elif t21 == n_triangle:
#            return t22
#        else:
#            return None
#    else:
#        temp = get_opposite_triangle(i, n_triangle - 6)
#        if temp == None:
#            return None
#        else:
#            return temp + 6



    
def make_integrand(dir1, dir2, triangle):
    return '(%s)*(%s)' % (grad_phi[triangle][dir1], grad_phi[triangle][dir2])
    #return ''.join(['(', grad_phi[triangle][dir1],')*(',grad_phi[triangle][dir2],')'])

def make_integrand_for_overlapped_tents(dir1, dir2, triangle1, triangle2):
    return '(%s)*(%s)' % (grad_phi[triangle1][dir1], grad_phi[triangle2][dir2])
    #return ''.join(['(', grad_phi[triangle1][dir1],')*(',grad_phi[triangle2][dir2],')'])

def get_full_integral(dir1, dir2, triangles):
    total = 0.0
    for triangle in triangles:
        expr = 'integrate(integrate(%s, %s), %s).evalf()' % (make_integrand(dir1, dir2, triangle), limits[triangle][0], limits[triangle][1])
        total += eval(expr)
    return np.float64(total)

def get_shared_triangles_integral(i, dir1, dir2, triangles):
    total = 0.0
    a = get_according_triangle(i, triangles)
    if len(a) == 0:
        return 0
    for t1, t2 in get_according_triangle(i, triangles):
        expr = 'integrate(integrate(%s, %s), %s).evalf()' % (make_integrand_for_overlapped_tents(dir1, dir2, t1, t2), limits[t1][0], limits[t1][1])
        total += eval(expr)
    return np.float64(total)
        

def detect_zone_quad(i):
    im, jm = divmod(i, m)
    i = im - 1
    j = jm - 1
    if (i < 0) or (j < 0) or (i >=n) or (j >= n):
        return None
    # full set of triangles
    triangles = set([])
    if i == 0:
        triangles |= set([3, 4, 5])
    elif i < n - 1:
        triangles |= set([0, 1, 2, 3, 4, 5])
    elif i == n - 1:
        triangles |= set([0, 1, 2])
    if j < 1:
        triangles -= set([2, 3, 4])
    elif j + 1 >= n:
        triangles  -= set([0, 1, 5])
    return frozenset(triangles)

def detect_zone_trapezium(i):
    im, jm = divmod(i, m)
    i = im - 1
    j = jm - 1
    if (i < 0) or (j < 0) or (i >= 2*n- 1) or (j >= n) or (i + j > 2*n-2):
        return None
    # full set of triangles
    triangles = set([])
    if i == 0:
        triangles |= set([3, 4, 5])
    elif i < n - 1:
        triangles |= set([0, 1, 2, 3, 4, 5])
    elif i == n - 1:
        triangles |= set([0, 1, 2, 3+6, 4+6, 5+6])
    else:
        triangles |= set([0+6, 1+6, 2+6, 3+6, 4+6, 5+6])
    if j < 1:
        triangles -= set([2, 3, 4, 2+6, 3+6, 4+6])
    elif j + 1 >= n:
        triangles  -= set([0, 1, 5, 0+6, 1+6, 5+6])
    if i + j == 2*n-2:
        triangles -= set([0+6, 4+6, 5+6])
    return frozenset(triangles)

def get_detector(sh_t):
    if sh_t == 1:
        return detect_zone_trapezium
    elif sh_t == 2:
        return detect_zone_quad
import sys
shape_type = int(sys.stdin.readline())

n = int(sys.stdin.readline())
mu1 = float(sys.stdin.readline())
e1 = float(sys.stdin.readline())
detect_zone = get_detector(shape_type)

g = e1/(2*(1+mu1))
m = n+2
k = 0
if shape_type == 1:
    k = m*(2*n+ 1)
elif shape_type == 2:
    k = m*m

f = open('test2.txt', 'w')
#print >>f, n, mu1, e1, k
print >>f, get_full_integral(0,1, frozenset([0,1,2,3,4,5]))
for i in range(6):
    print >>f, get_shared_triangles_integral(i, 0, 1, frozenset([0,1,2,3,4,5]))
f.close()    
#for i in xrange(k):
#    print >>f, i, detect_zone_quad(i)
#coeff1 = e1/(1-mu1**2)
#coeff2 = g / 2

#out = sys.stdout
coeffs_list = []
coeffs_map = [] # for points
zone_coeffs = dict() 
def get_coeffs():
    for i in xrange(k):
        zone = detect_zone(i)
        if zone == None :
            coeffs_map.append(None)
        elif not zone_coeffs.has_key(zone):
            b = [[[],[]],[[],[]]]
            b[0][0] = [e1/(1-mu1**2) * get_full_integral(0, 0, zone) + g/2 * get_full_integral(1, 1, zone)]
            b[0][0].extend([e1/(1-mu1**2) * get_shared_triangles_integral(i, 0, 0, zone) + g/2 * get_shared_triangles_integral(i, 1, 1, zone) for i in xrange(6)])
    	    
            b[0][1] = [e1/(1-mu1**2) * mu1 * get_full_integral(0, 1, zone) + g/2 * get_full_integral(1, 0, zone)]
            b[0][1].extend([e1/(1-mu1**2) * mu1 * get_shared_triangles_integral(i, 0, 1, zone) + g/2 * get_shared_triangles_integral(i, 1, 0, zone) for i in xrange(6)])
			
            b[1][0] = [e1/(1-mu1**2) * mu1 * get_full_integral(1, 0, zone) + g/2 * get_full_integral(0, 1, zone)]
            b[1][0].extend([e1/(1-mu1**2) * mu1 * get_shared_triangles_integral(i, 1, 0, zone) + g/2 * get_shared_triangles_integral(i,0, 1, zone) for i in xrange(6)])
            #b[1][0] = b[0][1]
	
            b[1][1] = [e1/(1-mu1**2) * get_full_integral(1, 1, zone) + g/2 * get_full_integral(0, 0, zone)]
            b[1][1].extend([e1/(1-mu1**2) * get_shared_triangles_integral(i, 1, 1, zone) + g/2 * get_shared_triangles_integral(i, 0, 0, zone) for i in xrange(6)])
            #b[1][1] = b[0][0]

            coeffs_list.append(b)
            ind = len(coeffs_list) - 1
            zone_coeffs[zone] = ind
            coeffs_map.append(ind)
        else:
            coeffs_map.append(zone_coeffs[zone])

get_coeffs()
#f.close()
print >>sys.stdout, len(coeffs_list)
for zone in xrange(len(coeffs_list)):
	for i in range (2):
		for j in xrange(2):
			for p in xrange(7):
				print >>sys.stdout, str(coeffs_list[zone][i][j][p])
for i in xrange(k):
	if coeffs_map[i] != None:
		print >>sys.stdout, str(coeffs_map[i])
	else:
		print >>sys.stdout, '-1'
print >>sys.stdout, 'Sucessful'

